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Abstract 

We study the joint effect of thermal bath fluctuations and an external noise tuning 
activity of cytotoxic cells on the triggered immune response in a growing cancerous 
tissue. The immune response is assumed to be primarily mediated by effector cells 
that develop a cytotoxic activity against the abnormal tissue. The kinetics of such 
a reaction is represented by an enzymatic-like Michaelis-Menten two step process. 
Effective free-energy surface for the process is further parameterised by the fluctu- 
ating energy barrier between the states of high and low concentration of cancerous 
cells. By analysing the far from equilibrium escape problem across the fluctuating 
potential barrier, we determine conditions of the most efficient decay kinetics of the 
cancer cell-population in the presence of dichotomously fluctuating concentration 
of cytotoxic cells. 
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1 Introduction 



Mathematical models of population growth based on nonlinear ordinary dif- 
ferential equations have been widely studied (see e.g. [1,2]), since despite their 
simplicity they can often capture the essence of complex biological interactions 
and explain characteristics of proliferation phenomena. However, biological 
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processes are not purely deterministic: systems existing in nature are subject 
to different types of natural noises. Generally, one can divide them into two 
types: (i)-internal noises (thermodynamic fluctuations) that display system 
parameters variability at the thermodynamic equilibrium, and (ii)-external 
(mostly conformational) noises, due solely to the effect of time-variability of 
the environment. In consequence, the information about dynamics of com- 
plex biological systems with underlying conformational transitions may be 
relevantly addressed by use of the mesoscopic approach in which the time- 
evolution of a collective coordinate is governed by stochastic differential equa- 
tions [3]. In such an approach, system kinetics is described in terms of an 
analogue to the Kramers scenario with the Brownian particle wandering over 
an effective free energy potential. In many cases, the potential experienced 
by the Brownian particle (or, correspondingly by a collective coordinate of 
the system) cannot be regarded as static but as affected by some random 
fluctuations with characteristic time scale that is comparable with one of the 
time scales governing the passage over the potential barrier separating basins 
of stationary states [4,5,6,7,8,9]. Some examples include the escape of O2 or 
CO ligand molecule from the myoglobin after dissociation, the ion channel 
kinetics in the lipid cell membrane or protein folding [6]. Another illustration 
of a coupling between the external noise and and a collective variable is an 
effect of cell-mediated immune surveillance against the cancer [10]. Most of 
tumoral cells bear antigens which are recognised as strange by the immune 
system. A response against these antigens may be mediated either by immune 
cells such as T-lymphocytes or other cells, not directly related to the immune 
system (like macrophages or natural killer cells). The process of damage to 
tumour proceeds via infiltration of the latter by the specialised cells which sub- 
sequently develop a cytotoxic activity against the cancer cell-population. The 
series of cytotoxic reactions between the cytotoxic cells and the tumour tissue 
may be considered to be well approximated [10] by a saturating, enzymatic-like 
process whose time evolution equations are similar to the standard Michaelis- 
Menten kinetics. The variability of kinetic parameters defining this process 
naturally affects the extinction of the tumour [10,11]. 

In this paper we investigate an extension of the above-mentioned problem: We 
discuss properties of a two state cancer-growth model subject to independent 
Markovian dichotomous noise (affecting the immunological response) and to 
an additive thermal noise (equilibrium fluctuations around stationary states 
representing low and high concentrations of cancer cells). As a model of the 
cancer cells dynamics, we have considered an overdamped Brownian particle 
moving in a two-well quasi-potential between reflecting and absorbing bound- 
aries, in the presence of noise that modulates the height of the barrier dividing 
two stable states of the population. Transitions from one state to the other 
(here: from a fixed-size tumour to a cancer-free state or vice versa ) are induced 
by additive thermal noise. In models of this kind, it is of particular interest 
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how fast the system escapes the potential well, surmounting the fluctuating 
potential barrier. In the case of the cancer growth model, we are especially 
interested in the rate at which the system escapes the fixed-size tumour state, 
that leads to the extinction of the tumour. The mean escape time from the 
meta-stable state across the fluctuating barrier may exhibit non-monotonic 
dependence on the characteristic time scale of these fluctuations, as has been 
exemplified in other studies [4,5,7,8,12,13,14,15]. Here, the phenomenon is dis- 
cussed for the model system in which existence of the dichotomous Markovian 
noise in one of the parameters may change relative stability of meta-stable 
states and therefore reverse the direction of the kinetic process. 
The following section presents briefly a generic model system that has been 
used for the analysis of cancer growth kinetics. Subsequent paragraphs are 
devoted to the presentation of the deterministic model and analysis of its 
stochastic counterpart. 



2 Generic Model System 



We adhere to the model of an overdamped Brownian particle moving in a 
potential field between absorbing and reflecting boundaries in the presence of 
noise which modulates the barrier height. The evolution of a state variable 
x(t) is described in terms of the Langevin equation 

dx dV(x) i _, . . . . . dV ± (x) i _, . 

* = ~ sr + "« f) + 9{xHt) = ~ & + "«*>■ (1) 

Here £(£) is a Gaussian process with zero mean and correlation 

< mm > = s(t - f) 



i.e. the Gaussian white noise of intensity a = y / 2k B T arising from the heat 
bath of temperature T), whereas r)(t) stands for a Markovian dichotomous 
noise switching between two levels {A + , A~} with mean frequency 7 and cor- 
relation time I/27. This means that its autocorrelation function is 

( (V(t) ~ (V)) (V(f) - (V)) ) = (A+ ~ A " )2 e -^l. 



Both noises are assumed to be statistically independent, i.e. < £(t)r)(s) >= 0. 
Eq. (1) can be considered as describing an overdamped motion of the state 
variable, subject to an effective force 

dV ± (x) dV(x) . 1 , , 
' +A ± g{x). 



dx dx 
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Based on eq. (1), one can write a set of the Fokker-Planck equations which 
describe the evolution of probability density of finding the state variable in a 
"position" x at time t: 



d t p(x,A ± ,t) = d x 



dV ± (x) 1 2o 



p(x, A ± ,t) 



dx 2 

- 1 p(x,A ± ,t) + 1P (x,A^,t) (2) 

In the above equations time has dimension of [length] 2 /energy due to a friction 
constant that has been "absorbed" in a time variable. 
With the initial condition 

p(V±,x s ,t)\ t=0 = ]-5(x-x s ), (3) 
the equations for the mean-first-passage time {MFPT) read: 



(4) 



1 +/ \ — / \ dV (x) dr (x) . a 2 d 2 r (x) 

where t + (x) and t~(x) denote MFPT for V + (x) and ^"(x), respectively. 

The overall MFPT for the system 

MFPT(x) = MFPT + (x) + MFPT~(x) = t + (x) + t~(x) (5) 

can be obtained after solving eq.(4) with appropriate boundary conditions. 
Here we assume a motion between a reflecting boundary a and an absorbing 
boundary b: 



dr+(x) 

dx 
dr~ (x) 
dx 



x=a 



x=a 








(6) 

r+(x)\ x=b = 

T~(x)\ x=b = 



According to a vast literature on the subject, the most interesting aspect of 
the escape over such a fluctuating potential is the non-monotonic behaviour of 
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the MFPT as a function of the driving noise [4,8,12,13,14]. In particular, we 
expect that for the frequency of potential switching 7 tending to zero (i.e. in 
the limit of a long correlation time of the dichotomous noise i](t)) the overall 
MFPT will be a mean value of MFPTs for both configurations, 

lim MFPT{V + ,V~,-f) = - [mFPT(V + ] + MFPT(V~)] , (7) 



whereas for 7 tending to infinity (in the limit of a short correlation time) the 
system will "experience" a mean barrier: 

Urn MFPT(V + ,V-,-f) = MFPT + ^-J , (8) 



with the MFPT(V + ) and MFPT(V~) obtained from formula 

MFPT( X ) = 1 /*,exp /.feexp (9) 



a; 



for V = l^, separately. 

Although the solution of (4) is usually unique, a closed, "ready to use" analyt- 
ical formula for r can be obtained only for the simplest cases of the potentials. 
More complex cases require either use of approximation schemes [7], pertur- 
bative approach [8] or direct numerical evaluation methods [9]. 



3 Population Model of Cancer Growth 



We use the predator-prey model ([3], [11]) to describe the cancer cells pop- 
ulation growth in presence of cytotoxic cells. The population dynamics can 
be described as follows: First, the cytotoxic cells bind to the tumour cells at 
a rate proportional to the kinetic constant k\\ second, the cancer cells which 
have been bound are killed and the complex dissociates at a rate proportional 
to k 2 . The process can be described schematically: 

X + Y Z Y + P . (10) 



Here X represents the population of tumour cells. Similarly, Y, Z and P rep- 
resent active cytotoxic cells, bound cells and dead tumour cells, respectively. 
In a given (small) volume element, there is an upper limit iV to the number of 
cells which may be present, given that each cell has a typical diameter equal 
to a. From now on, we will use normalised cellular densities: (x — jj instead 
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of X, etc). Following the original presentation [11], we assume that (i) cancer 
cells undergo replication at a rate proportional to the time constant A; (ii) 
as a result of cellular replication in limited volume, a diffusive propagation of 
cancer cells is possible, with transport coefficient Aa 2 ; (iii) dead cancer cells 
undergo elimination at a rate proportional to a certain constant (iv) local 
cytotoxic cell population remains constant, i.e. Y + Z = const; (v) free cy- 
totoxic cells can move with a "diffusion" coefficient D. The spatio-temporal 
evolution of the tumour due to the above processes can be then described by 
the set of balance equations: 

' § = A[l - (x + p)]x - hEyx + Aa 2 (l - p)V 2 x + Xa 2 xV 2 p 

|K = -kiyx + k 2 z + DV 2 y 
< at w 2 y (n) 

§f = hyx - k 2 z 
§ = k 2 Ez - hp 



where 

y + z = E = const. (12) 



In the limit when the effector cells diffuse much faster than the cancer cells 
propagate by cellular replication and in which the dead cells are rapidly elim- 
inated, the spatial distribution of Y and Z cells equilibrates rapidly with 
respect to the local density of living tumour cells and the above scheme of 
kinetics can be recasted in the form of the scalar problem [10,11]: 

?l = { l-ex)x-P^- + V 2 x ) (13) 
at x + 1 



where 



hx k 2 kiE 
x = — , 9=—, (3 = — — , t = Xt. 
k 2 ki A 



In our paper, we will consider the spatially homogeneous form of eq. (13). 
As presented here, the model has a long history of analytical studies ([11], 
[16]). In a slightly modified form it turned out to be also of practical use in 
biophysical modelling of radiation-induced damage production and processing 
[17,18,19]. In particular, the kinetic scheme eq.(10) has been adapted for the 
purpose of studying kinetics of double-strand breaks rejoining and formation 
of simple chromosome exchange aberrations [17] after DNA exposure to ion- 
ising radiation. Similar kinetics has been proposed in analysis of saturable 
repair models [18,19] devoted to study evolution of radio-biological damage. 
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In the latter, the saturable repair modelled by a Michaelis-Menten kinetics 
describes processing of damage by enzyme systems that can be overloaded. 
In the forthcoming section, we will analyse some versions of the model eq.(13) 
taking into account time variability of parameters as an effect of the environ- 
mental noise. Below, following [11], we briefly remind its deterministic prop- 
erties. 

Equation (13) can be considered as describing an overdamped motion of the 
state variable moving in a quasi - "free energy potential" : 

x 9t^ 

V(x) = -— + — + /3x-/3hi(x + l). (14) 

where V(x) has at most three extrema (stationary points of the system): 

xi=0, (15) 



1-9 + ^(1 + 9^-4(39 
X2 = - Q (16) 



and 



1-9 



x 3 



4(39 



29 



(17) 



Stability analysis reveals a strong dependence on 9: 

(i) For 9 > 1 and f3 > 1, Xi is the only minimum of V(x). For 9 > 1 and 
< (3 < 1, X\ becomes maximum and Xi is a new minimum. 

(ii) For 9 < 1, outside the region with < (3 < ^p-, the properties of the 
system are identical to (i). Inside the region, V(x) has two minima: x±, x 2 and 
one maximum at x 3 . For certain values of parameters, namely for 

9<1, 0</5<^£, (18) 



the system is bistable 1 . In a forthcoming section, we will consider this system 
subject to the joint effect of independent noises: to a multiplicative dichoto- 
mous noise in f3 with exponential time-correlation and to a white additive 
noise representing thermal fluctuations. 

1 Depending on 9, there exists a unique value of f3 = (3q, at which both states x\ 
and X3 are equally stable. Further stochastic analysis of the model is performed for 
9 = 0.1, and consequently (3q w 2.669 
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Fig. 1. Cancer growth kinetics has been studied for four indicated cases. MFPT was 
computed for intervals [a, b], where a and b are reflecting and absorbing boundaries, 
respectively. Arrows point directions of a transition. 

4 Stochastic Model of Cancer Growth 



We investigate the system described by the equation 

(i 1* IT 

_ = (1 - te)l - (/3 + A±)— + ^), 



(19) 



where is a two-state, Markovian noise and is the Gaussian noise of 
intensity a. In this form the model includes the influence of the heat bath 
(modelled by a memoryless Gaussian noise) and the fluctuations in the im- 
munological response of the organism (here assumed to be represented by a 
symmetric dichotomous noise in parameter (3). Correspondingly, the process 
of population growth and decay can be described as a motion of a fictitious 
particle in a potential switching between two conformational states (cf. Fig. 

!)• 



V ± (x) = -- + — + ([3±A)(x- ln(x + 1)). 



(20) 



For negligible additive noise and small concentration of cancerous cells, this 
model resembles standard Verhulst equation with a perturbing multiplica- 
tive dichotomous noise and as such has been shown [3] to exhibit a com- 
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2 4 6 8 10 



Fig. 2. Linear approximation of the original potential. The line connects minima of 
the potential V(x, A = 0) with the maximum of V(x, A^). Dotted lines: V(x, A^), 
solid lines: the approximated potential. 

plex scenario of noise-induced transitions observable in a pattern of the sta- 
tionary probability density. Here, we will address kinetic properties of this 
model by studying the mean first passage time (4) between high and low 
population states in the system. In order to compute MFPT, as a working 
example we use an approximation of the potential by a linear slope with 
a reflecting barrier placed at x\ = a = and and absorbing barrier at 
x 2 = b = [1 - 6 + ^(1 + 9) 2 - 4/30] (20) -1 . In a linear setup, it is frequently 
possible to formulate analytical solutions to the problem [5,12]. Moreover, as 
expected based on former analysis [15], due to the smoothing of the kinetics by 
irregular diffusive motions, the solutions for the piecewise linear potential are 
not qualitatively different from those for smooth differentiable potentials. In 
such fully analytic expression for the MFPT can be obtained, even 

though the algebra involved in such an evaluation requires use of symbolic 
computer softwares (here: Maple 7 procedure dsolve/numeric/BVP). 
Based on the original potential, we build its linear approximation in the fol- 
lowing way: In general, inclusion of a noise term in the parameter f3 influences 
positions of stationary states of the potential V(x, A 1 * 1 ) and affects their rel- 
ative stability. However, if the effective change in the barrier height is very 
small compared to its total height, the location of extrema can be considered 
constant (cf. Fig. (2)). Our analysis demonstrates that even in such cases of 
weak perturbation of the barrier height, presence of noise has dramatic con- 
sequences on the overall kinetics of cancer and its extinction. 

In order to keep constant the positions x±, x 2 and x 3 of the extrema when 
changing the barrier height V(x, A^), we approximate them with positions of 
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the extrema of the potential V(x,A = 0). Thus, x\ = 0, x 2 = 1 ^"^fl"^' ~ 

and x 3 = 1 6 ^^e^ ^ • Next, we connect points [xi, V(xi, A = 0)], [x 2 , V(x 2 , A + )] 
and [x 3 ,V(x 3 , A = 0)] or, respectively, [x 1 ,V(x 1 , A = 0)], [x 2 ,V(x 2 , A~)] and 
[x 3 , V(x 3 , A = 0)] with straight lines, as has been depicted in Fig.2. 
By defining 



A A := yl + 9 2 + 29 (1 — 2 (f3 + A)) (21) 

and 



A := A A=0 = ^1 + 6^ + 26(1-2(3), (22) 



the height of the barrier seen by the particle located close to the £3 can be 
expressed as 

V( I3 (A),A) - y (l2 (A),A) = Zd^!_^±i±ii) + 



Similarly, the height of the barrier calculated from the bottom of the left-side 
minimum x\ = is 



Piecewise linear approximation of the potential results in: 



and 



V( Z ,(0),A) = J 9 - 1 + W- 1 + ^ + 
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T p IT,*T lst IT p$t 




Fig. 3. Relative MFPTs as a function of the rate of the barrier fluctuations 7 
computed for various values of A, for an approximated linear potential. Parameters 
of the generic model has been set to (3 = 2.669 = (3q, a = 0.5 and 9 = 0.1. 

(26) 




For the purpose of analysis, the MFPT has been computed for 4 chosen cases 
(cf. Figure 1): (i) transition from large to small population when the small 
population state is a global minimum; (ii) transition from large to small pop- 
ulation when the large population state is a global minimum; (iii) transition 
from small to large population when the small population state is a global 
minimum, (iv) transition from small to large population when the large pop- 
ulation state is a global minimum. 

Figure 3 displays the graphs of MFPT as a function of the barrier fluctua- 
tion rate 7 and the noise intensity A. The results are presented in the form 
of a ratio T p /Ti of forward/backward mean first transition times calculated 
as a fraction of a corresponding ratio of transition times estimated for static 
(A = 0) barriers. 



The results indicate that the escape kinetics (from right to left and vice versa) 
depends on the slope of the approximated linear potential (cf. eqs. (23-26) and 
exhibits characteristic features of the noise-induced enhancement of the acti- 
vation process. Moreover, for the critical (deterministic) value of parameter 
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Fig. 4. Relative escape time as a function of the rate of the barrier fluctuations 
7 estimated for various duration of the integration time t t as compared with the 
down-hill relaxation time t r for a static barrier. 



(3 = (3o, fluctuations of the barrier may facilitate the forward/backward trans- 
fer and within given intervals of frequencies either the growth or diminishment 
of cancer population can be expected. For either of the transitions (from small- 
to large- and from large- to small concentrations of cancerous cells), depending 
on the frequency of the barrier fluctuations 7, there exists an optimal value of 
the barrier fluctuation rate, for which the mean first passage time is minimal 
[4,5,8,14,15]). The analysis reveals that the intensity of the effect depends on 
the strength of the barrier noise (the minimum/maximum of the relative T v jT\ 
deep ens/ lowers for larger A). Also, the resonant frequency of the MFPT min- 
imum 7 m j„ increases monotonically with the intensity a of the additive white 
Gaussian noise as displayed in Figure 6. 

Numerically estimated MFPT values have been compared with their asymp- 
totes for 7 — > 00, 7 — > evaluated directly from formulas (8), where 

x — b a 2 2A 2A 
MFPT(x) = — + ^[exp —(6 - a) - exp — (x - a)] (27) 



and 



V(a) - V ± {b) 
a — b 



(28) 
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Fig. 5. Same as in Figure 4, reverse transition over a fluctuating barrier, 
for MFPT(V ± ), or 

V(a) - V(b) 



A = 



a — b 



(29) 



for MFPT( V + v ), where a — x 1 and b = x 2 . The results display full agree- 
ment with numerical evaluation. 

Finally, we have used a recently proposed [8] approximated method for study- 
ing activation over a fluctuating barrier that involves considering separately 
the slow and fast components of barrier fluctuations and applies for any value 
of the correlation time r = I/7. By decomposing the barrier noise into two 
independent terms 



V (t)= Vs (t)+ Vf (t) (30) 

the escape problem may be considered as a three-dimensional Markovian pro- 
cess described by a joint probability distribution P(x,r] s ,r)f,t) which evolves 
accordingly to the Fokker-Planck equation with two separate operators re- 
sponsible for the time evolution of slow and fast variables. In particular, by 
definition, r) s (t) remains constant while the trajectory climbs the barrier and 
its dynamics can be analysed by the kinetic approach [8,12]. On the con- 
trary, rjf vanishes for the barrier correlation times slightly greater than the 
relaxation time from the top of the barrier to the bottom of the well, so 
that the probability distribution P(x,r) s ,r)f,t) can be deconvoluted to the 
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Fig. 6. Position of the MFPT minimum ^ m in as a function of noise intensity a 

form p(x,r]f,t;r] s )g(r] s ,t) and in consequence, the rate theory formalism can 
be safely applied. The equilibration process in a fast r]f(t) variable leads then 
to an effective Fokker-Planck equation (with a new form of the quasi-potential 
[8,13]) from which the appropriate effective MFPT(i] s ) can be calculated. For 
dichotomous switching the procedure captures ideology of the kinetic rate 
estimated as an inverse of 



UFPT ™ = 2[fc£ + +(*+ + + T hl - < 31 > 

where kinetic rates k + , k~ describe escape kinetics in two different configura- 
tions of the barrier. Note, that although the above formula resembles a typical 
kinetic scenario [12] which is known not to produce the resonant activation 
effect, the evaluation of the MFPT(r] s ) takes into account fast kinetics of the 
dynamics hidden in the form of kinetic constants evaluated from the effec- 
tive Fokker-Planck equation. Crucial for the approach is the determination of 
the value of the integration interval t t which is describing climbing stage of 
the process. Its duration for the unperturbed potential equals the relaxation 
time t r from the top to the bottom of the well. Obviously, fluctuations of the 
potential break that equality [12,8]. However, since we are not discussing the 
relationship between the process of climbing up and relaxing down the fluctu- 
ating barrier, but rather need a time-estimate for the following escape event, 
as a first approximation for t t the MFPT from the top of the barrier to the 
bottom of the well has been used. 



To test the method, we have presented in Figures 4 and 5 the exact analytical 
results for the triangular barrier (cf. Figure 2) along with the approximation 
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based on the approach described above. The agreement with the original re- 
sults is quite good and depends only slightly on the time of integration t t of 
the slow r)(t) component, as discussed elsewhere [8]. Nevertheless, we expect 
that such a comparison will be less positive for kinetics over a relatively low 
barriers (data not shown). 



5 Conclusions 

We have considered a standard problem of an escape over a fluctuating poten- 
tial barrier for a system describing efficiency of the immunological response 
against cancer. The basic component of the model is the mode of action of 
immunologically active cells that behave as catalysts in a chemical reaction. 
Analytical results for the model have been obtained after approximating a 
smooth double well potential by a piecewise-linear one. Although the method 
presented in this paper allows to estimate only the shortest possible time of 
transition (by placing an absorbing boundary condition at the top of the bar- 
rier), the results demonstrate fine sensitivity of the process to the correlation 
time of the barrier noise. In particular, by controlling the frequency of barrier 
fluctuations (or equivalently, the correlation time of changes in the effector 
cells response), the process of tissue growth can be reversed. Further analy- 
sis should tackle the problem of obvious nonlinearity of the process and the 
interplay of the noise with the stability criteria for the double-well system. 
The formalism of a partial noise-averaging method [8] may provide a way to 
analyse this case. 

It seems worthy mentioning that despite its simplicity, the model eq.(ll) can 
be considered as a generic one for kinetic description of many processes that 
use allosteric transition and phosphorylation of proteins in their metabolic 
pathways [17,18]. However, to ensure modelling of relevant interactions and to 
determine proper estimates of the kinetic parameters, an additional computer 
analysis of rate constants is required, which in turn has to rely on experimen- 
tal observations. In particular, local stability analysis of such models leads to 
conditions on rate parameters and steady state concentrations that need to be 
compared with experimental data. If the binding kinetics was controlable in 
such cases by an external (stochastic or deterministic) time-dependent field, 
the analysis similar to presented in this paper could be used for prediction of 
a most effective action of the catalyst leading to an expected shortest reaction 
time. In modelling tumor treatment or experiments on cell killing by external 
agents based on ordinary differential equations, a very common assumption 
is that treatment modifies the growth, so that the Verhulst- or Gompertz- 
type growth rate has to be accompanied by an extra term supressing the 
proliferating population [19]. The inhibition and deceleration of the tumor 
growth can often be assumed to follow a saturable repair model [2,17,18,19] of 
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the Michaelis-Menten kinetics. Therefore, consequences of the above analysis 
may appear realistic in those situations where, given a full identification of the 
parameters, the cellular populations would exhibit extinction controlled by ex- 
ternal, time dependent fields like time- dependent administration of cytotoxic 
agents or therapeutic radiation. 
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